1 Setup

1.2 Helper functions

# function for printing out html or latex tables 
print_table = function(data, format = "html", digits = 2){
  if(format == "html"){
    data %>% 
      kable(digits = digits) %>% 
      kable_styling()
  }else if(format == "latex"){
    data %>% 
      xtable(digits = digits) %>%
      print(include.rownames = F,
            booktabs = T)
  }
}

# root mean squared error
rmse = function(x, y){
  return(sqrt(mean((x - y)^2)))
}

actual_counterfactual_plot = function(x, ylabel) {
  p = ggplot(data = x, 
              mapping = aes(x = clipindex,
                            y = value,
                            fill = colorindex,
                            group = model)) +
    stat_summary(geom = "bar",
                 fun = "mean",
                 color = "black",
                 position = position_dodge(0.7),
                 width = 0.7) +
    stat_summary(geom = "linerange",
                 fun.data = "mean_cl_boot",
                 position = position_dodge(0.7),
                 size = 1) +
    geom_point(data = x %>% mutate(value = ifelse(model == "rating", value, NA)),
               position = position_jitterdodge(dodge.width = 0.7, jitter.width = 0.2),
               shape = 21,
               color = "black",
               size = 1,
               alpha = 0.2) +
    scale_fill_manual(values = c("red2", "green2", "black")) +
    facet_grid(index_actual ~ index_counterfactual) +
    geom_text(data = df.text %>%
                drop_na(label),
              mapping = aes(x = x, y = y, label = as.character(label)),
              size = 6,
              position = position_dodge(width = .7),
              hjust = 0.5) +
    coord_cartesian(xlim = c(0.5, 2.5),
                    ylim = c(-5, 105),
                    expand = F,
                    clip = "off") +
    labs(y = ylabel) +
    theme(text = element_text(size = 20),
          legend.position = "none",
          axis.text.x = element_blank(),
          axis.ticks.x = element_blank(),
          axis.title.y = element_text(margin = margin(0, 10, 0, 0)),
          panel.spacing.y = unit(0.8, "cm"),
          plot.margin = unit(c(0.2, 0.2, .8, .2), "cm"),
          axis.title.x = element_blank(),
          panel.grid = element_blank(),
          strip.background = element_blank(),
          panel.border = element_rect(colour = "black", fill = NA))
  # print(p)
}

actual_counterfactual_threeballs_plot = function(x, ylabel) {
  p = ggplot(data = x, 
              mapping = aes(x = ball,
                            y = value,
                            group = model,
                            fill = colorindex)) +
    stat_summary(fun = "mean",
                 geom = "bar",
                 color = "black",
                 position = position_dodge(0.9),
                 width = 0.9) +
    stat_summary(fun.data = "mean_cl_boot",
                 geom = "linerange",
                 color = "black",
                 position = position_dodge(0.9)) +
    geom_point(data = x %>% 
                 mutate(value = ifelse(model == "rating", value, NA)),
               position = position_jitterdodge(jitter.width = 0.3,
                                               jitter.height = 0,
                                               dodge.width = 0.9),
               shape = 21,
               color = "black",
               size = 1,
               alpha = 0.15) +
    facet_wrap(~clip, ncol = 8) +
    geom_text(data = df.text, 
              mapping = aes(x = x,
                            y = y,
                            label = as.character(label)),
              size = 5, 
              position = position_dodge(width = .9), 
              hjust = 0.5, 
              fontface = "bold") +
    coord_cartesian(xlim = c(0.5, 2.5), 
                    ylim = c(-5, 120), 
                    expand = F,
                    clip = "off") +
    scale_y_continuous(breaks = seq(0, 100, 25)) +
    scale_fill_manual(values = c("red2", "green2", "black")) +
    scale_color_manual(values = c("red2", "green2", "black")) +
    labs(y = ylabel) +
    theme(text = element_text(size = 16),
          legend.position = "none",
          axis.title.x = element_blank(),
          panel.spacing.y = unit(c(.3), "cm"),
          panel.grid = element_blank(),
          strip.background = element_blank(),
          strip.text = element_blank(),
          panel.border = element_rect(colour = "black", fill = NA))
  # print(p)
}

2 Experiment 1: 2 balls

2.1 Clips

2.2 Read in data

# causal judgments
df.exp1.causal = read.csv("../../data/experiment1_causal.csv", 
                          header = T,
                          stringsAsFactors = F)

# counterfactual judgments
df.exp1.counterfactual = read.csv("../../data/experiment1_counterfactual.csv",
                                  header = T,
                                  stringsAsFactors = F)

# Information about each clip 
df.exp1.clipinfo = tibble(clip = 1:18,
                          outcome_actual = c(0, 0, 0, 0, 0, 0, 0, 1, 0, 
                                             1, 1, 0, 1, 1, 1, 1, 1, 1),
                          outcome_counterfactual = c(0, 0, 0, 1, 1, 1, 0, 0, 1, 
                                                     0, 1, 1, 0, 0, 1, 0, 1, 1),
                          index_actual = rep(c("actual miss", 
                                               "actual close", 
                                               "actual hit"), each = 6),
                          index_counterfactual = rep(rep(c("counterfactual miss",
                                                           "counterfactual close", 
                                                           "counterfactual hit"),
                                                         each = 2),
                                                     3)) %>% 
  mutate(index_actual = factor(index_actual, levels = c("actual miss", 
                                                        "actual close", 
                                                        "actual hit")),
         index_counterfactual = factor(index_counterfactual, 
                                       levels = c("counterfactual miss",
                                                  "counterfactual close", 
                                                  "counterfactual hit")))

# Structure data 
df.exp1.causal.long = df.exp1.causal$data %>%
  enframe() %>% 
  mutate(participant = 1:n()) %>% 
  separate_rows(value, sep = ",") %>% 
  mutate(name = rep(c("clip", "rating"), n()/2),
         order = rep(rep(1:18, each = 2), max(participant)),
         value = as.numeric(value)) %>% 
  pivot_wider(names_from = name,
              values_from = value) %>% 
  arrange(participant, clip, rating) %>% 
  left_join(df.exp1.clipinfo,
            by = "clip")

df.exp1.counterfactual.long = df.exp1.counterfactual$data %>%
  enframe() %>% 
  mutate(participant = 1:n()) %>% 
  separate_rows(value, sep = ",") %>% 
  mutate(name = rep(c("clip", "rating"), n()/2),
         order = rep(rep(1:18, each = 2), max(participant)),
         value = as.numeric(value)) %>% 
  pivot_wider(names_from = name,
              values_from = value) %>% 
  arrange(participant, clip, rating) %>% 
  left_join(df.exp1.clipinfo,
            by = "clip")

2.3 Read in model predictions

2.4 Counterfactual condition

2.5 Causal condition

2.5.2 Stats

2.5.3 Tables

term estimate std.error lower 95% HDI upper 95% HDI
intercept -18.33 3.36 -23.81 -12.78
index_actualactualclose 4.31 3.40 -1.40 9.83
index_actualactualhit 4.06 4.84 -4.01 11.88
index_counterfactualcounterfactualclose -17.52 3.37 -22.96 -12.01
index_counterfactualcounterfactualhit -61.37 4.38 -68.54 -54.30
outcome_actual 92.36 5.11 84.04 100.96

3 Experiment 2: Brick & Teleport

3.1 Clips

3.2 Read in data

df.exp2.causal_first = read.csv("../../data/experiment2_causal_first.csv",
                                header = T,
                                stringsAsFactors = F)
df.exp2.counterfactual_first = read.csv("../../data/experiment2_counterfactual_first.csv",
                                        header = T,
                                        stringsAsFactors = F)

# Information about each clip

df.exp2.clipinfo = tibble(clip = 1:20,
                          outcome_actual = c(0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 
                                             0, 1, 1, 1, 1, 1, 1, 1, 0, 1),
                          outcome_counterfactual = c(0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 
                                                     1, 1, 0, 0, 1, 0, 1, 1, 1, 1),
                          index_actual = c(rep(c("actual miss", 
                                                 "actual close",
                                                 "actual hit"),
                                               each = 6),
                                           rep("practice", 2)),
                          index_counterfactual = c(rep(rep(c("counterfactual miss",
                                                             "counterfactual close",
                                                             "counterfactual hit"),
                                                           each = 2), 
                                                       3),
                                                   rep("practice", 2))) %>% 
  mutate(index_actual = factor(index_actual, levels = c("actual miss", 
                                                        "actual close", 
                                                        "actual hit")),
         index_counterfactual = factor(index_counterfactual, 
                                       levels = c("counterfactual miss",
                                                  "counterfactual close", 
                                                  "counterfactual hit")))

# Structure data 
df.exp2.causal_first.long = df.exp2.causal_first$data %>% 
  enframe() %>% 
  mutate(participant = 1:n()) %>% 
  separate_rows(value, sep = ",") %>% 
  mutate(name = rep(c("clip", "rating"), n()/2),
         order = rep(rep(1:40, each = 2), max(participant)),
         value = as.numeric(value)) %>% 
  pivot_wider(names_from = name,
              values_from = value) %>% 
  mutate(question = rep(rep(c("causal", "counterfactual"), each = 20), max(participant)),
         condition = "causal_first") %>% 
  left_join(df.exp2.clipinfo,
            by = "clip") %>%
  arrange(participant, question, clip) %>%
  select(participant, question, clip, rating, everything())

df.exp2.counterfactual_first.long = df.exp2.counterfactual_first$data %>% 
  enframe() %>% 
  mutate(participant = 1:n()) %>% 
  separate_rows(value, sep = ",") %>% 
  mutate(name = rep(c("clip", "rating"), n()/2),
         order = rep(rep(1:40, each = 2), max(participant)),
         value = as.numeric(value)) %>% 
  pivot_wider(names_from = name,
              values_from = value) %>% 
  mutate(question = rep(rep(c("counterfactual", "causal"), each = 20), max(participant)),
         condition = "counterfactual_first") %>% 
  left_join(df.exp2.clipinfo,
            by = "clip") %>%
  arrange(participant, question, clip) %>%
  select(participant, question, clip, rating, everything())

# combine data from both conditions
df.exp2.combined.long = df.exp2.causal_first.long %>%
  bind_rows(df.exp2.counterfactual_first.long %>% 
              mutate(participant = participant + 
                       max(df.exp2.causal_first.long$participant)))

3.3 Read in model predictions

3.4 Counterfactual condition

3.4.3 Tables

term estimate std.error lower 95% HDI upper 95% HDI
intercept 13.38 2.58 9.12 17.55
conditioncounterfactual_first -1.11 3.65 -6.99 4.82
index_actualactualclose -1.48 2.78 -6.02 3.11
index_actualactualhit -4.42 2.40 -8.32 -0.46
index_counterfactualcounterfactualclose 37.94 3.33 32.39 43.35
index_counterfactualcounterfactualhit 73.80 4.17 67.00 80.57
conditioncounterfactual_first:index_actualactualclose -0.55 3.84 -7.11 5.73
conditioncounterfactual_first:index_actualactualhit -0.85 3.41 -6.52 4.85
conditioncounterfactual_first:index_counterfactualcounterfactualclose 1.23 4.79 -6.62 9.18
conditioncounterfactual_first:index_counterfactualcounterfactualhit -0.02 5.92 -9.86 9.55

3.5 Causal condition

3.5.1 Plots

set.seed(1)

func_exp2_causal_plot = function(condition.name, model.name){

df.plot = df.exp2.combined.long %>%
  filter(question == "causal", clip <= 18) %>%
  filter(condition %in% condition.name) %>%
  mutate(rating = abs(rating),
         clipindex = rep(c(1, 2), nrow(.) / 2)) %>%
  left_join(df.exp2.causal.model,
            by = c("clip", 
                   "outcome_actual", 
                   "outcome_counterfactual", 
                   "index_actual", 
                   "index_counterfactual")) %>%
  pivot_longer(cols = c("rating", all_of(model.name)),
               names_to = "model",
               values_to = "value") %>% 
  mutate(model = factor(model, levels = c("rating", model.name))) %>%
  mutate(colorindex = ifelse(outcome_actual == 1, 2, 1),
         colorindex = ifelse(model != "rating", 3, colorindex),
         colorindex = as.factor(colorindex),
         index_actual = factor(index_actual, 
                               levels = c("actual miss", 
                                          "actual close", 
                                          "actual hit")),
         index_counterfactual = factor(index_counterfactual, 
                                       levels = c("counterfactual miss", 
                                                  "counterfactual close", 
                                                  "counterfactual hit")))

df.text = df.plot %>%
  expand(index_actual, index_counterfactual, clipindex, model) %>%
  mutate(label = rep(c("1 (A)", "2 (B)", "3", "4", "5 (A)", "6 (B)",
                       "7", "8 (C)", "9", "10", "11", "12 (C)",
                       "13 (D)", "14 (E)", "15", "16", "17 (D)", "18 (E)"), each = 2),
         label = ifelse(model != "rating", NA, label),
         y = -15,
         x = clipindex,
         colorindex = NA)

actual_counterfactual_plot(df.plot, "causal judgment")
}

plot.list = list(func_exp2_causal_plot(condition.name = "counterfactual_first",
                                       model.name = "combined"),
                 func_exp2_causal_plot(condition.name = "causal_first",
                                       model.name = "combined"))

# model.name = "combined"
# condition.name = "counterfactual_first"
# condition.name = "causal_first"
# 
# func_exp2_causal_plot(condition.name = condition.name,
#                       model.name = model.name),
# 
# ggsave(str_c("../../figures/plots/exp2_", condition.name ,"_causal_bars.pdf"),
#        width = 8,
#        height = 6)

3.5.2 Stats

3.5.3 Tables

3.5.3.1 Bayesian regression

term estimate std.error lower 95% HDI upper 95% HDI
intercept -24.07 3.50 -29.80 -18.27
conditioncounterfactual_first 11.47 4.91 3.29 19.40
index_actualactualclose 8.17 3.44 2.39 13.62
index_actualactualhit 13.42 4.90 5.42 21.30
index_counterfactualcounterfactualclose -30.20 3.72 -36.41 -24.36
index_counterfactualcounterfactualhit -50.29 4.60 -57.94 -42.84
outcome_actual 98.53 5.40 89.53 107.57
conditioncounterfactual_first:index_actualactualclose -5.39 4.94 -13.38 2.80
conditioncounterfactual_first:index_actualactualhit -16.97 7.02 -28.61 -5.75
conditioncounterfactual_first:index_counterfactualcounterfactualclose -8.13 5.17 -16.51 0.35
conditioncounterfactual_first:index_counterfactualcounterfactualhit -25.51 6.44 -36.53 -15.18
conditioncounterfactual_first:outcome_actual 10.72 7.60 -1.63 23.34

4 Experiment 3: 3 balls

4.1 Clips

4.2 Read in data

# clipinfo
df.exp3.clipinfo = tibble(clip = 1:32,
                          outcome_both = c(0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 
                                           0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1),
                          outcome_a = c(0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 
                                        0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1),
                          outcome_b = c(0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 
                                        0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1),
                          outcome_none = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                                           1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1),
                          clipindex = rep(1:2, 16))

# COUNTERFACTUAL JUDGMENTS 
df.exp3.counterfactual = read.csv("../../data/experiment3_counterfactual.csv",
                                  header = T,
                                  stringsAsFactors = F)

# demographics
df.exp3.counterfactual.demographics = df.exp3.counterfactual$extra %>% 
  enframe() %>% 
  separate_rows(value) %>% 
  mutate(value = as.numeric(value),
         name = rep(c("gender", "age", "difficulty", "smoothness", "time", "condition"),
                    n()/6),
         participant = rep(1:(n()/6), each = 6)) %>% 
  pivot_wider(names_from = name,
              values_from = value) %>% 
  mutate(condition = factor(condition, levels = 20:21, labels = c("A", "B")),
         gender = factor(gender, levels = 1:2, labels = c("female", "male"))) %>% 
  select(participant, condition, gender, age, time, smoothness, difficulty)

# judgments
df.exp3.counterfactual.long = df.exp3.counterfactual$data %>%
  enframe() %>% 
  separate_rows(value) %>% 
  mutate(value = as.numeric(value),
         name = rep(c("clip", "rating"), n()/2),
         participant = rep(1:(n()/64), each = 64),
         order = rep(rep(1:32, each = 2), n()/64)) %>% 
  pivot_wider(names_from = name,
              values_from = value) %>% 
  left_join(df.exp3.counterfactual.demographics %>% 
              select(participant, ball = condition),
            by = "participant") %>% 
  select(participant, ball, clip, everything()) %>% 
  arrange(participant, clip)

# CAUSAL JUDGMENTS
df.exp3.causal = read.csv("../../data/experiment3_causal.csv",
                          header = T,
                          stringsAsFactors = F)

# demographics 
df.exp3.causal.demographics = df.exp3.causal$extra %>% 
  enframe() %>% 
  separate_rows(value) %>% 
  mutate(value = as.numeric(value),
         name = rep(c("gender", "age", "difficulty", "smoothness", "time", 
                      "counterbalance", "replayType", "experiment"),
                    n()/8),
         participant = rep(1:(n()/8), each = 8)) %>% 
  pivot_wider(names_from = name,
              values_from = value) %>% 
  mutate(counterbalance = ifelse(participant %in% c(1, 3, 5, 6, 9), 1, counterbalance),
         gender = factor(gender, levels = 1:2, labels = c("female", "male"))) %>% 
  select(participant, gender, age, counterbalance, time, smoothness, difficulty)

# judgments
df.exp3.causal.long = df.exp3.causal$data %>%
  enframe() %>% 
  separate_rows(value) %>% 
  mutate(value = as.numeric(value),
         name = rep(c("clip", "x", "y", "z", "rating1", "rating2"), n()/6),
         participant = rep(1:(n()/(32*6)), each = (32*6)),
         order = rep(rep(1:32, each = 6), n()/(32*6))) %>% 
  filter(!name %in% c("x", "y", "z")) %>% 
pivot_wider(names_from = name,
            values_from = value) %>% 
  left_join(df.exp3.causal.demographics %>% 
              select(participant, counterbalance),
            by = "participant") %>% 
  mutate(A = ifelse(counterbalance == 1, rating1, rating2),
         B = ifelse(counterbalance == 1, rating2, rating1)) %>% 
  select(-rating1, -rating2) %>% 
  pivot_longer(cols = c(A, B),
               names_to = "ball",
               values_to = "rating") %>% 
  select(participant, clip, ball, rating, order) %>% 
  arrange(participant, clip, ball)

4.3 Read in model predictions

4.3.3 Heuristic model

l.features = fromJSON("data/features.json")

df.features.initial = l.features[["appearance"]] %>%
  cbind(l.features[["initialVelocity"]][["ballA"]]) %>%
  cbind(l.features[["initialVelocity"]][["ballB"]]) %>%
  cbind(l.features[["initialVelocity"]][["ballE"]]) %>%
  setNames(c(str_c("ball", LETTERS[c(1, 2, 5)], "_appearance"),
             "ballA_velx",
             "ballA_vely",
             "ballB_velx",
             "ballB_vely",
             "ballE_velx",
             "ballE_vely")) %>%
  mutate(clip = 1:nrow(.)) %>%
  pivot_longer(cols = starts_with("ball")) %>%
  separate(name, into = c("ball", "property")) %>%
  pivot_wider(names_from = property, values_from = value) %>% 
  mutate(ball = str_remove_all(ball, pattern = "ball"))

# information about collisions
df.features.collisions = data.frame()
ballnames = c("ballA", "ballB", "ballE")

for (i in 1:32) {
  df.tmp = data.frame()
  for (j in 1:length(ballnames)) {
    ncollisions = length(l.features[["collisions"]][[ballnames[j]]][["object"]][[i]])
    if (ncollisions > 0) {
      tmp = data.frame(
        ball = ballnames[j],
        object = l.features[["collisions"]][[ballnames[j]]][["object"]][[i]],
        time = l.features[["collisions"]][[ballnames[j]]][["time"]][[i]],
        pre.velx = l.features[["collisions"]][[ballnames[j]]][["pre"]][[i]][["x"]],
        pre.vely = l.features[["collisions"]][[ballnames[j]]][["pre"]][[i]][["y"]],
        post.velx = l.features[["collisions"]][[ballnames[j]]][["post"]][[i]][["x"]],
        post.vely = l.features[["collisions"]][[ballnames[j]]][["post"]][[i]][["y"]],
        ncollision = 1:ncollisions
      )
      df.tmp = df.tmp %>%
        rbind(tmp)
    }
  }
  df.features.collisions = df.features.collisions %>%
    rbind(df.tmp %>%
            mutate(clip = i) %>%
            select(clip, ball, everything()))
}

# find end of clip
tmp = df.features.collisions %>%
  group_by(clip) %>%
  filter(ball == "ballE",
         str_detect(object, "goal|Left|Right")) %>%
  group_by(clip) %>%
  mutate(endclip = max(time)) %>%
  select(clip, endclip) %>%
  ungroup() %>%
  distinct()

# remove events that happen after the end of the clip
df.features.collisions = df.features.collisions %>%
  left_join(tmp, by = "clip") %>%
  mutate(endclip = ifelse(is.na(endclip), 400, endclip)) %>%
  group_by(clip) %>%
  filter(time <= endclip)

# E initially at rest or moving?
tmp.movingE = df.features.initial %>%
  filter(ball == "E") %>%
  mutate(E.moving = ifelse(velx == 0 & vely == 0, 1, 0)) %>%
  select(clip, E.moving)

# contact with ball E
tmp.contactE = df.features.collisions %>%
  filter(object == "ballE") %>%
  mutate(ball = factor(ball, levels = c("ballA", "ballB"), labels = c("A", "B"))) %>%
  mutate(contactE = 1) %>%
  select(clip, ball, contactE) %>%
  group_by(clip, ball) %>%
  summarize(contactE = any(contactE %>% as.logical()) * 1) %>%
  left_join(expand.grid(clip = 1:32, ball = c("A", "B")), .,
            by = c("clip", "ball")) %>%
  mutate(contactE = ifelse(is.na(contactE), 0, contactE))

# number of collisions
tmp.ncollisions = df.features.collisions %>%
  filter(str_detect(object, "ball"),
         ball != "ballE") %>%
  mutate(ball = factor(ball, levels = c("ballA", "ballB"), labels = c("A", "B"))) %>%
  group_by(clip, ball) %>%
  count() %>%
  rename(ncollision = n)

# exclusivity (no contact with ball E by the other ball)
tmp.exclusive = df.features.collisions %>%
  filter(str_detect(object, "ball"),
         ball != "ballE") %>%
  count(clip, ball, object) %>%
  mutate(AE = ifelse(ball == "ballA" & object == "ballE", 1, 0),
         BE = ifelse(ball == "ballB" & object == "ballE", 1, 0)) %>%
  group_by(clip) %>%
  summarize(AE = any(AE %>% as.logical()) * 1,
            BE = any(BE %>% as.logical()) * 1) %>%
  mutate(A = ifelse(AE == 1 & BE == 0, 1, 0),
         B = ifelse(AE == 0 & BE == 1, 1, 0)) %>%
  select(clip, A, B) %>%
  pivot_longer(cols = -clip,
               names_to = "ball",
               values_to = "exclusive")

# impact
func_angle = function(x, y) {
  dot.prod = x %*% y
  norm.x = norm(x, type = "2")
  norm.y = norm(y, type = "2")
  theta = acos(dot.prod / (norm.x * norm.y))
  as.numeric(theta)
}

# impact on ball E 
tmp.impactE = df.features.collisions %>%
  rowwise() %>%
  filter(str_detect(object, "ball"),
         ball == "ballE") %>%
  mutate(pre.speed = abs(pre.velx) + abs(pre.vely),
         post.speed = abs(post.velx) + abs(post.vely),
         speed.diff = post.speed - pre.speed,
         direction.diff = func_angle(c(pre.velx, pre.vely), c(post.velx, post.vely)),
         direction.diff = ifelse(is.na(direction.diff),
                                 abs(atan2(post.vely - pre.vely, post.velx - pre.velx)),
                                 direction.diff)) %>%
  ungroup() %>%
  group_by(clip, object) %>%
  summarize(speed.diff = sum(speed.diff),
            direction.diff = sum(direction.diff)) %>%
  rename(ball = object) %>%
  mutate(ball = factor(ball, levels = c("ballA", "ballB"), labels = c("A", "B"))) %>%
  left_join(expand.grid(clip = 1:32, ball = c("A", "B")), .,
            by = c("clip", "ball")) %>%
  mutate_at(.vars = vars(speed.diff, direction.diff),
            .funs = ~ ifelse(is.na(.), 0, .)) %>%
  arrange(clip, ball) %>%
  rename(E.speed.diff = speed.diff,
         E.direction.diff = direction.diff)

# total impact 
tmp.impactTotal = df.features.collisions %>%
  rowwise() %>%
  filter(str_detect(object, "ballA|ballB")) %>%
  mutate(pre.speed = abs(pre.velx) + abs(pre.vely),
         post.speed = abs(post.velx) + abs(post.vely),
         speed.diff = post.speed - pre.speed,
         direction.diff = func_angle(c(pre.velx, pre.vely), c(post.velx, post.vely)),
         direction.diff = ifelse(is.na(direction.diff),
                                 abs(atan2(post.vely - pre.vely, post.velx - pre.velx)),
                                 direction.diff)) %>%
  ungroup() %>%
  group_by(clip, object) %>%
  summarize(speed.diff = sum(speed.diff),
            direction.diff = sum(direction.diff)) %>%
  rename(ball = object) %>%
  mutate(ball = factor(ball, levels = c("ballA", "ballB"), labels = c("A", "B"))) %>%
  left_join(expand.grid(clip = 1:32, ball = c("A", "B")), .,
            by = c("clip", "ball")) %>%
  mutate_at(.vars = vars(speed.diff, direction.diff),
            .funs = ~ ifelse(is.na(.), 0, .)) %>%
  arrange(clip, ball) %>%
  rename(total.speed.diff = speed.diff,
         total.direction.diff = direction.diff)

# transfer of force 
tmp.transfer = df.features.collisions %>%
  filter(str_detect(ball, "ballA|ballB"),
         str_detect(object, "ball")) %>%
  mutate(AE = ifelse(ball == "ballA" & object == "ballE", 1, NA),
         BE = ifelse(ball == "ballB" & object == "ballE", 1, NA),
         AB = ifelse((ball == "ballA" & object == "ballB"), 1, NA)) %>%
  mutate_at(.vars = vars(AE, BE, AB), .funs = ~ . * time) %>%
  filter(!is.na(AE) | !is.na(BE) | !is.na(AB)) %>%
  arrange(clip, time) %>%
  group_by(clip) %>%
  summarize(A = any(!is.na(AE)),
            B = any(!is.na(BE)),
            A = ifelse(any(!is.na(AB)) &
                         any(!is.na(BE)) &
                         max(AB, na.rm = T) < min(BE, na.rm = T),
                       T,
                       A),
            B = ifelse(any(!is.na(AB)) &
                         any(!is.na(AE)) &
                         max(AB, na.rm = T) < min(AE, na.rm = T),
                       T,
                       B)) %>%
  pivot_longer(cols = -clip,
               names_to = "ball",
               values_to = "transfer") %>% 
  arrange(clip, ball)

# collect features
df.features = df.features.initial %>%
  filter(ball != "E") %>%
  mutate(ball = factor(ball)) %>%
  mutate(moving = ifelse(velx == 0 & vely == 0, 0, 1),
         speed = abs(velx) + abs(vely)) %>%
  left_join(tmp.contactE) %>%
  left_join(tmp.impactE) %>%
  left_join(tmp.impactTotal) %>%
  left_join(tmp.transfer %>% mutate(transfer = transfer * 1)) %>%
  left_join(tmp.movingE) %>%
  left_join(tmp.exclusive) %>%
  select(-c(appearance, velx, vely))

# remove temporary variables
rm(list = ls()[str_detect(ls(), "tmp")])

4.4 Counterfactual condition

4.5 Causal condition

4.5.1 Stats

4.5.1.2 Model comparison

No problematic observations found. Returning the original 'loo' object.
No problematic observations found. Returning the original 'loo' object.
No problematic observations found. Returning the original 'loo' object.
                        elpd_diff se_diff
fit.brm_exp3_causal_whs    0.0       0.0 
fit.brm_exp3_causal_wh  -125.2      15.6 
fit.brm_exp3_causal_w   -426.2      31.1 

4.5.1.3 Heuristic model

  • z-scored predictors and invidiual responses
 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: rating ~ moving + speed + contact_e + e_speed_diff + e_direction_diff + total_speed_diff + total_direction_diff + transfer + e_moving + exclusive + (1 | participant) 
   Data: df.regression.features %>% select(-c(clip, ball, o (Number of observations: 2624) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Group-Level Effects: 
~participant (Number of levels: 41) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     8.53      1.22     6.44    11.18 1.00     1163     1936

Population-Level Effects: 
                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept               49.79      1.46    46.91    52.74 1.00     1129
moving                   0.22      0.22     0.01     0.81 1.00     3399
speed                    2.06      0.85     0.44     3.67 1.00     2319
contact_e                0.38      0.35     0.01     1.29 1.00     3201
e_speed_diff             0.12      0.12     0.00     0.46 1.00     4435
e_direction_diff         1.04      0.74     0.04     2.71 1.00     2202
total_speed_diff         2.21      0.95     0.41     4.17 1.00     2208
total_direction_diff     3.99      0.91     2.14     5.68 1.00     2718
transfer                15.59      0.79    14.06    17.14 1.00     3555
e_moving                 0.18      0.17     0.01     0.62 1.00     3296
exclusive                4.39      0.73     2.97     5.81 1.00     3526
                     Tail_ESS
Intercept                1854
moving                   1547
speed                    1266
contact_e                2020
e_speed_diff             2268
e_direction_diff         1685
total_speed_diff         1393
total_direction_diff     1742
transfer                 2992
e_moving                 1738
exclusive                2538

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    33.21      0.46    32.33    34.10 1.00     5313     2845

Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

4.5.1.6 Individual participant regressions

df.fit = df.exp3.causal.long %>% 
  left_join(df.exp3.model.causal,
            by = c("clip", "ball"))

# priors 
priors = c(set_prior("normal(0,10)", class = "b", lb = 0))

# initial model fits (used for compilation)
fit.brm_exp3_causal_individual_baseline = 
  brm(formula = rating ~ 1,
      data = df.fit %>% 
        filter(participant == 1),
      save_all_pars = T,
      cores = 2,
      seed = 1,
      file = str_c("cache/brm_exp3_causal_individual_baseline"))

fit.brm_exp3_causal_individual_w = 
  brm(formula = rating ~ 1 + whether,
      data = df.fit %>% 
        filter(participant == 1),
      prior = priors,
      save_all_pars = T,
      cores = 2,
      seed = 1,
      file = str_c("cache/brm_exp3_causal_individual_w"))

fit.brm_exp3_causal_individual_wh = 
  brm(formula = rating ~ 1 + whether + how,
      data = df.fit %>% 
        filter(participant == 1),
      prior = priors,
      save_all_pars = T,
      cores = 2,
      seed = 1,
      file = str_c("cache/brm_exp3_causal_individual_wh"))

fit.brm_exp3_causal_individual_whs = 
  brm(formula = rating ~ 1 + whether + how + sufficient,
      data = df.fit %>% 
        filter(participant == 1),
      prior = priors,
      save_all_pars = T,
      cores = 2,
      seed = 1,
      file = str_c("cache/brm_exp3_causal_individual_whs"))

# update model fits for different participants
df.exp3.causal.individual_fit = df.fit %>%
  group_by(participant) %>%
  nest() %>%
  ungroup() %>% 
  # fit model for each participant
  mutate(fit_baseline = map(.x = data,
                            .f = ~ update(fit.brm_exp3_causal_individual_baseline,
                                          newdata = .x)),
         fit_w = map(.x = data,
                     .f = ~ update(fit.brm_exp3_causal_individual_w,
                                   newdata = .x)),
         fit_wh = map(.x = data,
                      .f = ~ update(fit.brm_exp3_causal_individual_wh,
                                    newdata = .x)),
         fit_whs = map(.x = data,
                       .f = ~ update(fit.brm_exp3_causal_individual_whs,
                                     newdata = .x))) %>% 
  mutate(fit_baseline = map(.x = fit_baseline,
                            ~ add_criterion(.x, criterion = c("loo", "waic"))),
         fit_w = map(.x = fit_w,
                     ~ add_criterion(.x, criterion = c("loo", "waic"))),
         fit_wh = map(.x = fit_wh,
                      ~ add_criterion(.x, criterion = c("loo", "waic"))),
         fit_whs = map(.x = fit_whs,
                       ~ add_criterion(.x, criterion = c("loo", "waic"))),
         r_w = map2_dbl(.x = data,
                        .y = fit_w,
                        .f = ~ cor(.x$rating, .y %>% 
                                     fitted() %>% 
                                     .[, 1])),
         r_wh = map2_dbl(.x = data,
                         .y = fit_wh,
                         .f = ~ cor(.x$rating, .y %>% 
                                      fitted() %>% 
                                      .[, 1])),
         r_whs = map2_dbl(.x = data,
                          .y = fit_whs,
                          .f = ~ cor(.x$rating, .y %>% 
                                       fitted() %>% 
                                       .[, 1])),
         rmse_w = map2_dbl(.x = data,
                        .y = fit_w,
                        .f = ~ rmse(.x$rating, .y %>% 
                                     fitted() %>% 
                                     .[, 1])),
         rmse_wh = map2_dbl(.x = data,
                         .y = fit_wh,
                         .f = ~ rmse(.x$rating, .y %>% 
                                      fitted() %>% 
                                      .[, 1])),
         rmse_whs = map2_dbl(.x = data,
                          .y = fit_whs,
                          .f = ~ rmse(.x$rating, .y %>% 
                                       fitted() %>% 
                                       .[, 1])),
         model_comparison = pmap(.l = list(baseline = fit_baseline,
                                           w = fit_w,
                                           wh = fit_wh,
                                           whs = fit_whs),
                                 .f = ~ loo_compare(..1, ..2, ..3, ..4)),
         best_model = map_chr(.x = model_comparison,
                              .f = ~ rownames(.) %>% .[1]),
         best_model = factor(best_model,
                             levels = c("..1", "..2", "..3", "..4"),
                             labels = c("baseline", "w", "wh", "whs")))

save(list = c("df.exp3.causal.individual_fit"),
     file = "data/exp3_causal_individual_fit.RData")

4.5.2 Plots

4.5.2.2 Selection

# \u2713 = check mark 
# \u2717 = cross mark 

check = "\u2713"
cross = "\u2717"

x_labels = c(str_c("\nW ", cross, "\nH ", check, "\nS ", cross),
             str_c("\nW ", check, "\nH ", check, "\nS ", check),
             str_c("\nW ", cross, "\nH ", cross, "\nS ", cross),
             str_c("\nW ", check, "\nH ", cross, "\nS ", cross),
             str_c("\nW ", check, "\nH ", check, "\nS ", cross),
             str_c("\nW ", check, "\nH ", check, "\nS ", cross),
             str_c("\nW ", cross, "\nH ", check, "\nS ", check),
             str_c("\nW ", cross, "\nH ", check, "\nS ", check),
             str_c("\nW ", cross, "\nH ", check, "\nS ", check),
             str_c("\nW ", cross, "\nH ", cross, "\nS ", cross))
  

df.plot = df.exp3.causal.long %>%
  filter(clip %in% c(3, 7, 15, 16, 23)) %>%
  group_by(clip, ball) %>%
  summarize(mean = mean(rating),
            low = smean.cl.boot(rating)[2],
            high = smean.cl.boot(rating)[3]) %>%
  left_join(df.exp3.causal.regression %>% select(clip, ball, w, wh, whs),
            by = c("clip", "ball")) %>%
  ungroup() %>%
  pivot_longer(cols = c(mean, w, wh, whs),
               names_to = "index",
               values_to = "value") %>% 
  mutate_at(.vars = vars(low, high), .funs = ~ ifelse(index == "mean", ., NA)) %>%
  mutate(index = factor(index, levels = c("mean", "w", "wh", "whs")),
         clip = factor(clip, levels = c(7, 23, 3, 15, 16),
                       labels = c("causal chain", "double prevention", "joint causation",
                                  "overdetermination", "preemption")))

df.labels = df.plot %>% 
  distinct(clip, ball) %>% 
  arrange(clip, ball) %>% 
  mutate(labels = x_labels,
         value = -10,
         index = NA)

func_load_image = function(clip){
  readPNG(str_c("../../figures/diagrams/exp3_clip", clip, ".png"))
}

df.clips = df.plot %>% 
  distinct(clip) %>% 
  arrange(clip) %>% 
  mutate(number = c(7, 23, 3, 15, 16),
         grob = map(.x = number, .f = ~ func_load_image(clip = .x)),
         index = NA,
         value = NA,
         ball = NA,
         label = str_c("clip ", number))

ball_a = readPNG("../../figures/diagrams/ballA.png")
ball_a = rasterGrob(ball_a, interpolate = TRUE)

ball_b = readPNG("../../figures/diagrams/ballB.png")
ball_b = rasterGrob(ball_b, interpolate = TRUE)

p = ggplot(data = df.plot, 
       mapping = aes(x = ball,
                     y = value,
                     group = index,
                     fill = index)) +
  geom_bar(stat = "identity", color = "black",
           position = position_dodge(0.9),
           width = 0.9) +
  geom_errorbar(mapping = aes(ymin = low, ymax = high),
                width = 0,
                size = 1,
                position = position_dodge(0.9)) +
  annotation_custom(grob = ball_a, xmin = 0.5, xmax = 1.5, ymin = -25, ymax = -7) + 
  annotation_custom(grob = ball_b, xmin = 1.5, xmax = 2.5, ymin = -25, ymax = -7) +
  geom_text(data = df.labels,
            mapping = aes(label = labels,
                          y = -Inf),
            vjust = 1.2,
            family = "Arial Unicode MS",
            size = 8) + 
  geom_custom(data = df.clips,
              mapping = aes(data = grob, x = 1.5, y = Inf),
              grob_fun = function(x) rasterGrob(x,
                                                interpolate = T,
                                                vjust = -0.25)) + 
  geom_text(data = df.clips,
            mapping = aes(label = label,
                          y = Inf,
                          x = -Inf),
            size = 9,
            hjust = -0.2,
            vjust = -4) + 
  facet_grid(~clip) +
  scale_fill_grey(start = 1,
                  end = 0,
                  labels = c("mean rating  ", expression(CSM[W] ~ " "),
                             expression(CSM[WH] ~ " "),
                             expression(CSM[WHS] ~ " "))) +
  labs(y = "causal responsibility", fill = "") +
  coord_cartesian(clip = "off") + 
  theme_bw() +
  theme(legend.text.align = 0,
        text = element_text(size = 20),
        panel.grid = element_blank(),
        legend.position = "bottom",
        axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        legend.text = element_text(size = 30),
        strip.text = element_text(size = 30),
        legend.spacing = unit(0.5, "cm"),
        legend.background = element_rect(fill = "transparent"),
        legend.margin = margin(t = 5, unit = "cm"),
        plot.margin = margin(t = 8, l = 0.2, r = 0.2, b = 0.1, unit = "cm"))
print(p)

# ggsave("../../figures/plots/exp3_selection_bars.pdf",
#        width = 20,
#        height = 10,
#        device = cairo_pdf)

4.5.2.3 Scatter plot

func_scatterplot = function(model){
  if(model == "heuristic"){
    xlabel = "Heuristic"
  }else{
    xlabel = bquote(CSM[.(toupper(model))])
  }
  
  l.models = list(w = fit.brm_exp3_causal_w, 
                  wh = fit.brm_exp3_causal_wh,
                  whs = fit.brm_exp3_causal_whs,
                  heuristic = fit.brm_exp3_causal_heuristic)
  
  df.data = df.exp3.causal.means %>% 
    left_join(df.exp3.model.causal, by = c("clip", "ball")) %>% 
    left_join(df.regression.features, by = c("clip", "ball"))
  
  df.plot = l.models[[model]] %>%
    fitted(newdata = df.data,
           re_formula = NA) %>% 
    as_tibble() %>% 
    clean_names() %>%
    bind_cols(df.data)
  
  p = ggplot(data = df.plot,
             mapping = aes(x = estimate, 
                           y = rating_mean)) +
    geom_abline(intercept = 0,
                slope = 1,
                linetype = 2) + 
    geom_smooth(aes(y = estimate,
                    ymin = q2_5,
                    ymax = q97_5),
                stat = "identity",
                color = "black") +
    geom_linerange(size = 0.75, 
                   mapping = aes(ymin = rating_low,
                                 ymax = rating_high),
                   color = "gray50") +
    geom_point(size = 2) + 
    scale_color_manual(values = c("black", "#e41a1c"),
                       guide = F) +
    coord_fixed(xlim = c(0, 100),
                ylim = c(0, 100)) +
    labs(x = xlabel,
         y = "responsibility judgment") +
    annotate(geom = "text",
             label = str_c(
               "r = ", cor(df.plot$estimate, df.plot$rating_mean) %>% 
                 round(2) %>% 
                 as.character() %>% 
                 str_sub(start = 2),
               "\nRMSE = ", rmse(df.plot$estimate, df.plot$rating_mean) %>% 
                 round(2)),
             x = 5,
             y = 95,
             size = 6,
             hjust = 0) +
    scale_x_continuous(breaks = seq(0, 100, 25)) +
    scale_y_continuous(breaks = seq(0, 100, 25)) +
    theme(text = element_text(size = 20))
  return(p)
}

# for creating and saving an individual scatter plot 
# model = "w"
# model = "wh"
# model = "whs"
# model = "heuristic"

# func_scatterplot(model)
# ggsave(str_c("../../figures/plots/exp3_", model, "_scatter.pdf"),
#        width = 5,
#        height = 5)

list(func_scatterplot(model = "w"),
     func_scatterplot(model = "wh"),
     func_scatterplot(model = "whs"),
     func_scatterplot(model = "heuristic")) %>%
  wrap_plots(ncol = 2)

4.5.2.4 Individual participant variance for a selection of clips (clustered)

# \u2713 = check mark 
# \u2717 = cross mark 

check = "\u2713"
cross = "\u2717"

x_labels = c(str_c("\nW ", cross, "\nH ", check, "\nS ", cross),
             str_c("\nW ", check, "\nH ", check, "\nS ", check),
             str_c("\nW ", cross, "\nH ", cross, "\nS ", cross),
             str_c("\nW ", check, "\nH ", cross, "\nS ", cross),
             str_c("\nW ", check, "\nH ", check, "\nS ", cross),
             str_c("\nW ", check, "\nH ", check, "\nS ", cross),
             str_c("\nW ", cross, "\nH ", check, "\nS ", check),
             str_c("\nW ", cross, "\nH ", check, "\nS ", check),
             str_c("\nW ", cross, "\nH ", check, "\nS ", check),
             str_c("\nW ", cross, "\nH ", cross, "\nS ", cross))

df.plot = df.exp3.causal.long %>%
  filter(clip %in% c(7, 23, 3, 15, 16)) %>%
  mutate(clip = factor(clip, levels = c(7, 23, 3, 15, 16),
    labels = c("causal chain", "double prevention", "joint causation",
      "overdetermination", "preemption")))

df.cluster = df.exp3.causal.individual_fit %>%
  select(participant, fit_whs) %>%
  mutate(estimates = map(fit_whs, tidy)) %>%
  select(participant, estimates) %>%
  unnest(estimates) %>%
  filter(str_detect(term, "b_"),
         !str_detect(term, "Intercept")) %>%
  mutate(term = str_remove(term, "b_")) %>%
  select(participant, term, estimate) %>%
  pivot_wider(names_from = term,
              values_from = estimate) %>% 
  mutate(group = ifelse(how > whether, "how", "whether"))

df.plot = df.plot %>% 
  left_join(df.cluster %>% 
              select(participant, group),
            by = "participant")

df.labels = df.plot %>% 
  distinct(clip, ball) %>% 
  arrange(clip, ball) %>% 
  mutate(labels = x_labels,
         rating = -10,
         group = NA,
         participant = NA)

func_load_image = function(clip){
  readPNG(str_c("../../figures/diagrams/exp3_clip", clip, ".png"))
}

df.clips = df.plot %>% 
  distinct(clip) %>% 
  arrange(clip) %>% 
  mutate(number = c(7, 23, 3, 15, 16),
         grob = map(.x = number, .f = ~ func_load_image(clip = .x)),
         group = NA,
         participant = NA,
         ball = NA,
         label = str_c("clip ", number))

ball_a = readPNG("../../figures/diagrams/ballA.png")
ball_a = rasterGrob(ball_a, interpolate = TRUE)

ball_b = readPNG("../../figures/diagrams/ballB.png")
ball_b = rasterGrob(ball_b, interpolate = TRUE)

# ggplot(df.plot,aes(x=ball,y=rating,group=id,color = best))+
p = ggplot(data = df.plot, 
       mapping = aes(x = ball,
                     y = rating,
                     group = participant,
                     shape = group)) +
  geom_line(mapping = aes(linetype = "individual",
                          color = group),
            size = 1,
            alpha = 0.3) +
  geom_point(mapping = aes(color = group),
             size = 1,
             alpha = 0.3) +
  stat_summary(mapping = aes(group = group,
                             color = group,
                             linetype = "mean"),
               fun = "mean",
               geom = "line",
               size = 1.5) +
  stat_summary(data = df.plot %>% filter(group == "how"),
               mapping = aes(group = group,
                             color = group),
               fun.data = "mean_cl_boot",
               geom = "pointrange",
               size = 1,
               shape = 19) +
  stat_summary(data = df.plot %>% filter(group == "whether"),
               mapping = aes(group = group,
                             color = group),
               fun.data = "mean_cl_boot",
               geom = "pointrange",
               size = 1,
               shape = 19) +
  annotation_custom(grob = ball_a, xmin = 0.5, xmax = 1.5, ymin = -30, ymax = -10) + 
  annotation_custom(grob = ball_b, xmin = 1.5, xmax = 2.5, ymin = -30, ymax = -10) +
  geom_text(data = df.labels,
            mapping = aes(label = labels,
                          y = -Inf),
            vjust = 1.2,
            family = "Arial Unicode MS",
            size = 8) + 
  geom_custom(data = df.clips,
              mapping = aes(data = grob, x = 1.5, y = Inf),
              grob_fun = function(x) rasterGrob(x,
                                                interpolate = T,
                                                vjust = -0.25)) +
  geom_text(data = df.clips,
            mapping = aes(label = label,
                          y = Inf,
                          x = -Inf),
            size = 9,
            hjust = -0.2,
            vjust = -4) + 
  facet_wrap(~clip,
             ncol = 8) +
  coord_cartesian(xlim = c(0.9, 2.1),
                  ylim = c(-5, 105),
                  expand = F,
                  clip = "off") +
  scale_y_continuous(breaks = seq(0, 100, 25),
                     labels = seq(0, 100, 25)) +
  scale_color_brewer(palette = "Set1",
                     guide = "none") +
  labs(y = "causal responsibility rating",
       linetype = "legend",
       shape = "") +
  theme_bw() +
  guides(linetype = guide_legend(override.aes = list(alpha = c(0.3, 1)),
                                 keywidth = unit(1.2, "cm")),
         shape = guide_legend(override.aes = list(color = c("#e41a1c", "#377eb8"),
                                                  shape = c(19, 19),
                                                  alpha = c(1, 1),
                                                  size = c(5, 5)))) +
  theme(panel.grid = element_blank(),
        text = element_text(size = 20),
        legend.position = c(0.505, 0.23),
        axis.title.x = element_blank(),
        legend.text = element_text(size = 20),
        legend.title = element_text(size = 24),
        legend.background = element_rect(fill = "transparent"),
        strip.text = element_text(size = 30),
        axis.text.x = element_blank(),
        legend.key = element_rect(fill = "transparent"),
        legend.box = "horizontal",
        legend.spacing.x = unit(0.1, "cm"),
        panel.spacing.x = unit(1, "cm"),
        plot.margin = margin(b = 5, l = 0.2, r = 0.2, t = 7.5, unit = "cm"))

print(p)

# ggsave(str_c("../../figures/plots/exp3_individual_variance_selection_lines_clustered.pdf"),
#        plot = p,
#        width = 20,
#        height = 8.5,
#        device = cairo_pdf)

4.5.3 Tables

4.5.3.2 Population level predictors in the mixed effects models

model term estimate std.error lower 95% HDI upper 95% HDI
CSM_w intercept 31.68 1.61 29.02 34.37
CSM_w whether 39.30 2.16 35.76 42.82
CSM_wh intercept 10.29 1.56 7.80 12.83
CSM_wh whether 28.40 2.68 23.96 32.70
CSM_wh how 36.07 2.60 31.80 40.27
CSM_whs intercept 10.43 1.59 7.82 13.01
CSM_whs whether 27.93 2.60 23.67 32.18
CSM_whs how 26.47 2.91 21.76 31.35
CSM_whs sufficient 30.90 3.06 25.77 35.96

4.5.3.3 Individual participants regression results

participant r_w r_wh r_whs rmse_w rmse_wh rmse_whs best_model
1 0.29 0.38 0.48 30.00 29.03 27.74 whs
2 0.20 0.77 0.77 39.32 27.94 27.55 whs
3 0.37 0.78 0.79 38.90 28.00 26.95 whs
4 0.41 0.55 0.63 43.62 40.60 38.31 whs
5 0.45 0.53 0.54 39.23 37.32 36.81 whs
6 0.18 0.54 0.54 40.84 36.32 36.07 whs
7 0.49 0.61 0.64 32.96 29.93 29.15 whs
8 0.39 0.63 0.68 38.37 33.26 31.64 whs
9 0.40 0.72 0.76 35.86 28.34 26.39 whs
10 0.28 0.52 0.53 29.13 26.28 25.85 whs
11 0.51 0.56 0.60 33.66 32.19 31.17 whs
12 0.56 0.57 0.70 32.78 32.08 28.84 whs
13 0.56 0.71 0.74 29.83 25.40 24.21 whs
14 0.43 0.47 0.50 33.72 32.82 32.23 whs
15 0.19 0.35 0.37 39.77 38.29 37.94 whs
16 0.62 0.67 0.76 40.49 37.98 34.39 whs
17 0.42 0.66 0.71 28.57 23.82 22.42 whs
18 0.31 0.35 0.39 37.78 37.09 36.57 whs
19 0.37 0.60 0.62 33.41 29.23 28.62 whs
20 0.57 0.61 0.68 36.47 35.00 32.71 whs
21 0.43 0.50 0.57 30.95 29.61 28.32 whs
22 0.51 0.67 0.71 33.72 29.36 28.00 whs
23 0.33 0.45 0.52 38.76 37.00 35.57 whs
24 0.46 0.63 0.69 34.01 29.96 28.22 whs
25 0.61 0.66 0.70 30.60 28.91 27.55 whs
26 0.57 0.71 0.76 40.43 35.03 32.80 whs
27 0.46 0.52 0.66 43.12 41.49 38.12 whs
28 0.30 0.45 0.46 29.33 27.62 27.38 whs
29 0.48 0.60 0.65 38.93 36.03 34.48 whs
30 0.18 0.70 0.70 33.06 24.95 24.71 whs
31 0.48 0.59 0.69 37.54 34.72 31.76 whs
32 0.37 0.61 0.70 40.62 35.94 32.92 whs
33 0.50 0.56 0.64 38.64 36.79 34.85 whs
34 0.24 0.51 0.52 44.50 40.86 40.19 whs
35 0.15 0.32 0.30 34.56 33.36 33.37 wh
36 0.19 0.84 0.82 45.38 29.58 29.70 wh
37 0.55 0.63 0.67 36.04 33.36 32.11 whs
38 0.52 0.61 0.67 36.03 33.39 31.53 whs
39 0.47 0.77 0.77 29.61 21.75 21.58 whs
40 0.46 0.78 0.79 32.41 23.60 22.69 whs
41 0.32 0.50 0.58 32.90 30.38 28.69 whs

4.5.3.4 CSMwhs predictions for selection of cases

clip ball difference whether how sufficient robust
7 A cmark (1) xmark (0.34) cmark (1) xmark (0) xmark (0.25)
7 B cmark (1) cmark (1) cmark (1) cmark (0.67) cmark (0.6)
23 A xmark (0.05) xmark (0) xmark (0) xmark (0) xmark (0)
23 B cmark (0.91) cmark (0.79) xmark (0) xmark (0) cmark (0.72)
16 A cmark (1) xmark (0.23) cmark (1) cmark (1) xmark (0.35)
16 B xmark (0) xmark (0) xmark (0) xmark (0) xmark (0)
3 A cmark (1) cmark (0.88) cmark (1) xmark (0.12) cmark (0.76)
3 B cmark (1) cmark (0.89) cmark (1) xmark (0.11) cmark (0.75)
15 A cmark (1) xmark (0.01) cmark (1) cmark (0.99) xmark (0.1)
15 B cmark (1) xmark (0.01) cmark (1) cmark (1) xmark (0.1)

4.5.3.5 CSMwhs predictions for all cases

clip ball outcome_both outcome_a outcome_b outcome_none difference whether how sufficient robust w wh whs heuristic rating
1 A 0 0 0 0 100 40 100 23 36 47 58 55 57 42
1 B 0 0 0 0 100 15 100 16 9 38 51 46 54 37
2 A 0 0 0 0 57 12 0 0 10 36 14 14 25 21
2 B 0 0 0 0 18 0 0 0 0 32 10 10 24 19
3 A 1 0 0 0 100 88 100 12 76 66 71 65 72 76
3 B 1 0 0 0 100 89 100 11 75 67 72 65 72 75
4 A 1 0 0 0 100 78 100 4 78 62 69 60 58 63
4 B 1 0 0 0 100 95 100 15 57 69 73 68 54 78
5 A 0 0 1 0 100 90 100 0 47 67 72 62 47 71
5 B 0 0 1 0 100 0 100 0 0 32 46 37 68 22
6 A 0 0 1 0 100 59 100 16 35 55 63 59 53 73
6 B 0 0 1 0 100 18 100 6 14 39 52 44 53 22
7 A 1 0 1 0 100 34 100 0 25 45 56 46 70 59
7 B 1 0 1 0 100 100 100 67 60 71 75 86 64 79
8 A 1 0 1 0 0 0 0 0 0 32 10 10 25 7
8 B 1 0 1 0 100 100 100 100 100 71 75 96 84 92
9 A 0 1 0 0 0 0 0 0 0 32 10 10 14 8
9 B 0 1 0 0 100 100 100 0 100 71 75 65 78 90
10 A 0 1 0 0 77 18 0 0 22 39 15 16 15 23
10 B 0 1 0 0 98 79 0 0 63 63 33 33 21 55
11 A 1 1 0 0 100 70 100 77 68 59 66 80 71 93
11 B 1 1 0 0 0 0 0 0 0 32 10 10 16 4
12 A 1 1 0 0 100 82 100 74 83 64 70 83 57 77
12 B 1 1 0 0 100 0 100 12 24 32 46 41 53 37
13 A 0 1 1 0 67 34 0 0 35 45 20 20 14 8
13 B 0 1 1 0 70 35 0 0 35 46 20 20 21 64
14 A 0 1 1 0 97 91 0 0 59 67 36 36 21 22
14 B 0 1 1 0 91 77 0 0 51 62 32 32 20 18
15 A 1 1 1 0 100 1 100 99 10 32 47 68 71 76
15 B 1 1 1 0 100 1 100 100 10 32 47 68 71 76
16 A 1 1 1 0 100 23 100 100 35 41 53 74 80 92
16 B 1 1 1 0 0 0 0 0 0 32 10 10 14 4
17 A 0 0 0 1 100 19 100 37 18 39 52 54 66 69
17 B 0 0 0 1 100 0 100 36 17 32 46 48 65 46
18 A 0 0 0 1 100 11 100 40 17 36 49 52 55 63
18 B 0 0 0 1 100 7 100 37 9 34 48 50 56 66
19 A 1 0 0 1 100 74 100 7 65 61 67 60 55 53
19 B 1 0 0 1 100 72 100 7 65 60 67 59 55 49
20 A 1 0 0 1 100 92 100 8 72 68 72 65 57 41
20 B 1 0 0 1 100 88 100 4 53 66 71 63 56 71
21 A 0 0 1 1 100 47 100 40 45 50 60 63 58 80
21 B 0 0 1 1 100 9 100 21 10 35 49 46 59 18
22 A 0 0 1 1 100 100 100 89 83 71 75 92 48 60
22 B 0 0 1 1 100 8 100 0 15 35 49 39 53 42
23 A 1 0 1 1 5 0 0 0 0 32 10 10 15 3
23 B 1 0 1 1 91 79 0 0 72 63 33 33 22 39
24 A 1 0 1 1 100 66 100 4 63 58 65 57 57 44
24 B 1 0 1 1 100 94 100 22 79 69 73 70 54 73
25 A 0 1 0 1 100 25 100 21 26 42 53 50 69 43
25 B 0 1 0 1 100 74 100 54 65 61 67 74 56 73
26 A 0 1 0 1 100 6 100 3 9 34 48 40 60 39
26 B 0 1 0 1 100 87 100 35 54 66 71 72 46 69
27 A 1 1 0 1 100 97 100 52 97 70 74 80 67 80
27 B 1 1 0 1 0 0 0 0 0 32 10 10 17 6
28 A 1 1 0 1 100 90 100 22 80 67 72 69 79 89
28 B 1 1 0 1 0 0 0 0 0 32 10 10 12 5
29 A 0 1 1 1 100 58 100 24 44 55 63 61 66 47
29 B 0 1 1 1 100 63 100 24 38 57 64 62 54 67
30 A 0 1 1 1 100 57 100 29 49 54 63 62 63 58
30 B 0 1 1 1 100 46 100 24 39 50 60 57 63 56
31 A 1 1 1 1 100 2 100 4 3 32 47 39 63 44
31 B 1 1 1 1 100 4 100 4 4 33 47 39 53 46
32 A 1 1 1 1 0 0 0 0 0 32 10 10 16 5
32 B 1 1 1 1 100 75 100 66 73 61 68 78 65 71

4.5.3.6 Heuristic model

term estimate std.error lower 95% HDI upper 95% HDI
intercept 49.79 1.46 47.39 52.23
moving 0.22 0.22 0.01 0.66
speed 2.06 0.85 0.64 3.43
contact_e 0.38 0.35 0.02 1.07
e_speed_diff 0.12 0.12 0.01 0.37
e_direction_diff 1.04 0.74 0.08 2.43
total_speed_diff 2.21 0.95 0.67 3.80
total_direction_diff 3.99 0.91 2.45 5.45
transfer 15.59 0.79 14.30 16.90
e_moving 0.18 0.17 0.01 0.51
exclusive 4.39 0.73 3.16 5.56

5 Session info

R version 3.6.2 (2019-12-12)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Mojave 10.14.6

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] forcats_0.5.0    stringr_1.4.0    dplyr_0.8.5      purrr_0.3.3     
 [5] readr_1.3.1      tidyr_1.0.2      tibble_2.1.3     tidyverse_1.3.0 
 [9] patchwork_1.0.0  egg_0.4.5        gridExtra_2.3    png_0.1-7       
[13] tidybayes_2.0.2  Hmisc_4.3-1      ggplot2_3.3.0    Formula_1.2-3   
[17] survival_3.1-11  lattice_0.20-40  janitor_1.2.1    jsonlite_1.6.1  
[21] xtable_1.8-4     corrr_0.4.2      brms_2.12.0      Rcpp_1.0.4      
[25] broom_0.5.5      lme4_1.1-21      Matrix_1.2-18    kableExtra_1.1.0
[29] knitr_1.28      

loaded via a namespace (and not attached):
  [1] readxl_1.3.1         backports_1.1.5      plyr_1.8.6          
  [4] igraph_1.2.5         splines_3.6.2        svUnit_0.7-12       
  [7] crosstalk_1.1.0.1    rstantools_2.0.0     inline_0.3.15       
 [10] digest_0.6.25        htmltools_0.4.0      rsconnect_0.8.16    
 [13] fansi_0.4.1          magrittr_1.5         checkmate_2.0.0     
 [16] cluster_2.1.0        modelr_0.1.6         matrixStats_0.56.0  
 [19] xts_0.12-0           prettyunits_1.1.1    jpeg_0.1-8.1        
 [22] colorspace_1.4-1     rvest_0.3.5          haven_2.2.0         
 [25] xfun_0.12            callr_3.4.2          crayon_1.3.4        
 [28] zoo_1.8-7            glue_1.3.2           gtable_0.3.0        
 [31] webshot_0.5.2        pkgbuild_1.0.6       rstan_2.19.3        
 [34] abind_1.4-5          scales_1.1.0         mvtnorm_1.1-0       
 [37] DBI_1.1.0            miniUI_0.1.1.1       viridisLite_0.3.0   
 [40] htmlTable_1.13.3     HDInterval_0.2.0     foreign_0.8-76      
 [43] stats4_3.6.2         StanHeaders_2.21.0-1 DT_0.13             
 [46] htmlwidgets_1.5.1    httr_1.4.1           threejs_0.3.3       
 [49] arrayhelpers_1.1-0   RColorBrewer_1.1-2   ellipsis_0.3.0      
 [52] acepack_1.4.1        farver_2.0.3         pkgconfig_2.0.3     
 [55] loo_2.2.0            nnet_7.3-13          dbplyr_1.4.2        
 [58] labeling_0.3         tidyselect_1.0.0     rlang_0.4.5         
 [61] reshape2_1.4.3       later_1.0.0          cellranger_1.1.0    
 [64] munsell_0.5.0        tools_3.6.2          cli_2.0.2           
 [67] generics_0.0.2       ggridges_0.5.2       evaluate_0.14       
 [70] fastmap_1.0.1        yaml_2.2.1           fs_1.3.2            
 [73] processx_3.4.2       nlme_3.1-145         mime_0.9            
 [76] xml2_1.2.5           compiler_3.6.2       bayesplot_1.7.1     
 [79] shinythemes_1.1.2    rstudioapi_0.11      reprex_0.3.0        
 [82] stringi_1.4.6        highr_0.8            ps_1.3.2            
 [85] Brobdingnag_1.2-6    nloptr_1.2.2.1       markdown_1.1        
 [88] shinyjs_1.1          vctrs_0.2.4          pillar_1.4.3        
 [91] lifecycle_0.2.0      bridgesampling_1.0-0 data.table_1.12.8   
 [94] httpuv_1.5.2         R6_2.4.1             latticeExtra_0.6-29 
 [97] bookdown_0.18        promises_1.1.0       boot_1.3-24         
[100] colourpicker_1.0     MASS_7.3-51.5        gtools_3.8.1        
[103] assertthat_0.2.1     withr_2.1.2          shinystan_2.5.0     
[106] parallel_3.6.2       hms_0.5.3            rpart_4.1-15        
[109] coda_0.19-3          minqa_1.2.4          snakecase_0.11.0    
[112] rmarkdown_2.1        shiny_1.4.0.2        lubridate_1.7.4     
[115] base64enc_0.1-3      dygraphs_1.1.1.6